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Abstract 

The galactic supershells are cavities in the interstellar medium. These shells can 
be explained by introducing the concept of superbubbles, the theoretical result of 
multiple supernova. The superbubbles can be analytically described if the ambient 
medium has a constant density both in the so-called bursting phase and in the sub- 
sequent adiabatic expansion. In order to solve the expansion of superbubbles in the 
ISM, which is a non-homogeneous medium, a numerical technique is used that di- 
vides the sphere into many sectors. By varying the time of the bursting phenomenon 
and the time over which the phenomenon is followed, elliptical and hour-glass shapes 
or vertical walls can be obtained. Application of the developed theory/code to the 
super-shell associated with GW 46.4-1-5 and with GSH 238 allows us to say that the 
suggested physical parameters are consistent with our theory. The map of the expand- 
ing superbubble's velocity can be tentatively traced by generating random points on 
the expanding surface. The structure of the galactic plane as a result of the evolution 
of many super-bubbles was simulated by adopting the percolation theory in order to 
generate new OB associations. 

Key words: ISM: bubbles - ISM: clouds - Galaxy: disk - galaxies: starburst - 
methods: numerical 



1. Introduction 

The term supershell was observationally defined by Heiles (1979) where eleven H I objects 
were examined. The supershells have been observed as expanding shells, or holes, in the HI- 
column density distribution of our galaxy, in the Magellanic Clouds, in the dwarf irregular 
galaxy Ho-II (Puche et al. 1992), and in the nearby dwarf galaxy IC 2574 (Walter et al. 
1998). The dimensions of these objects span from 100 pc to 1700 pc and often present elliptical 
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shapes or elongated features, which are difficult to explain based on an expansion in a uniform 
medium. These structures are commonly explained through introducing theoretical objects 
named bubble or superbubble; these are created by mechanical energy input from stars (see 
for example Pikel'ner 1968 ; Weaver et al. 1977). Thus, the origin of a supershell is not 
necessarily a superbubble, and other possible origins include collisions of high- velocity clouds 
(see for example Tenorio-Tagle, Bodenheimer 1988 ; Santillan et al. 1999 ). The worm is 
another observed feature that may, or may not, be associated with a wall of a superbubble. 
Galactic worms were ffist identified as irregular, vertical columns of atomic gas stretching from 
the galactic plane; now, similar structures are found in radio continuum and infrared maps (see 
for example Koo et al. (1992)). 

The models that explain supershells as being due to the combined explosions of supernova 
in a cluster of massive stars are now briefly reviewed. Semi-analytical and hydrodynamical 
calculations are generally adopted. In semi-analytical calculations the thin-shell approximation 
can be the key to obtaining the expansion of the superbubble; see, for example, McCray 1987; 
McCray , Kafatos 1987; Mac Low, McCray 1988; Igumenshchev et al. 1990; Basu et al. 1999. 

The thin-shell approximation has already been used in a variety of different problems, 
and with both analytical and numerical approaches (aside from the review by Bisnovatyi-Kogan, 
Silich 1995, see Begelman, Li 1992; Moreno et al. 1999). Thus the validity and limitations 
of the method are well known. For instance, modelling is fast and simple because only shell 
dynamics are included in an approximate way, while fluid variables are not included at all. 
The price that one has to pay is a lack of knowledge of the density and velocity profiles and, 
obviously, the onset and development of turbulence, instabilities or mixing cannot be followed. 
These facts limit the applications of this method to derive only general shapes, approximate 
expansion rates, and gross features of the surface mass-density distributions of the shells. A 
detailed comparison with the results of so-called thin-layer approximation, as described for 
example in Silich (1992) is now being carried out. 

• The two main assumptions, which are (i) all swept-up gas accumulates infinitely in a thin 
shell just after the shock front and (ii) the gas pressure is uniform inside the cavity, are in 
common with the two models. 

• The numerical equations of motion are different. For examples (1) and (5) are for during 
the bursting phase, and Silich's set from (5) to (9) is for the three Cartesian coordinates. 

• In our model the shear due to galactic rotation is neglected in a ffist approximation. See 
subsubsection 2.7.1 for a late addition on shear. 

• Our treatment includes the concept of efficiency of the simulation, equation (54), which is 
neglected by Silich. 

Another source of comparison can be made with the 3D models of Silich et al. (1996), in which 
a number of processes not considered here are now included: 



2 



• Two mechanisms of mass injection into the cavity are adopted: thermal evaporation of 
both the cold expanding shell, and the clouds engulfed by the shell, and dynamic disruption 
of clouds penetrating into the bubble interior. 

• A power-law density and temperature distribution inside the superbubble during the entire 
evolution. 

• A cooling function that takes into account the metallicity gradients in the Galaxy. 

However the astrophysical results of Silich et al. (1996) are basically similar to ours, 
and the maximum elongation of the superbubbles along the z-direction can reach 2 kpc in 
both models. Similar results are also obtained concerning the morphologies: in both models 
the hourglass- shape is obtained when equahty between the bursting time and the elapsed time 
is assumed. The vertical walls considered here are the belt of subsection 4.2 of Silich et al. 
(1996). 

The hydrodynamical approximation was used by MacLow et al. (1989). As for the effect 
of magnetic fields, a semi-analytical method was introduced by (Ferriere et al. 1991) and a 
magneto-hydrodynamic code has been adopted by various authors (Tomisaka 1992; Tomisaka 
1998; Kamaya 1998). 

The plane of the numerical analysis is now outlined. 

A numerical code that solves the momentum equation coupled with the variation of 
pressure in the presence of the injection of mechanical luminosity and adiabatic losses was 
developed in section 2; testing with numerical hydro-dynamic calculation was developed in 
subsubsection 3.1. 

The analytical "expansion law" of the superbubbles can be set up both in the bursting 
phase ( see section 2) and in the adiabatic phase after the SN burst stops, (see subsection 2.1). 

The various types of obtained shape (elliptical,egg,V, vertical and hourglass) are de- 
scribed in subsection 3.3. 

A numerical law that describes the raising of the vertical walls is deduced in subsubsec- 
tion 3.3.1. 

The map of the expanding superbubble's velocity can be tentatively traced by generating 
random points on the expansion surface ( see subsection 4.1). 

The developed theory was then applied to the supershells associated with 
GW 46.4+5 and GSH-238 (see subsection 4.1 and subsection 4.2, respectively). 

The H I structure in the galactic plane can be tentatively simulated by using the perco- 
lation framework to trace the spiral-arm structure ( see section 5 and appendix 2). 

By adopting the framework of radiation coming from the layer comprised between two 
shells, we can deduce the center darkening law (see appendix 1). 
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2. The Approximation Used 



In our case, the starting equation for the evolution of the superbubble (McCray 
1987; McCray , Kafatos 1987) is momentum conservation applied to a pyramidal section, 
characterised by a solid angle, AQj-. 

|(AM,i?,)=pi?|A^, (1) 

where the pressure of the surrounding medium is assumed to be negligible and the mass is 
confined into a thin shell with mass AMj. The subscript j was added here in order to note 
that this is not a spherically symmetric system. Due to the fact that p is uniform in the cavity 
of the superbubble, a summation to obtain the total volume, such as 

V^ = ^Afi,i?5/3, (2) 
is necessary to determine its value. The mass conservation equation for a thin shell is 

AM, = IrPjAQj. (3) 
The pressure is enclosed in the energy conservation equation, 

' (4) 



7 — 1 c/t dt 

where L = £'o-RsN/4vr denotes the mechanical luminosity injected into a unit solid angle and 
7 =5/3. Equation (4) can be expanded, obtaining 
dp _ iv(7- 1) p dV 

dt V dt' ^ ^ 

Formulae (1) and (5) will be our basic equations to numerically integrate in the bursting phase. 
An approximation concerning the pressure can be obtained by ignoring the second term of the 
rhs in (4); this leads to 

where t is the considered time, i?sN the rate of supernova explosions, Eq the energy of each 
supernova, and V is computed as in equation (2). On continuing to consider the case of constant 
density, the volume becomes: 

V = fR\ (7) 
Equations (1) and (6) lead to 

dt^P^^)-^^R- 
R^R = AR°' is imposed to integrate this equation. After adopting the initial condition of i? = 
at t = and assuming p is constant irrespective of R or t, the equation on the expansion speed 
is obtained 
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^^ a + 13EoRs. ,_ (9) 



By integrating the previous equation, "the expansion law" is obtained, 



R 



Ana 



1/5 / p p \ 1/5 



P 



(10) 



which is identical to equation (10.34) of McCray (1987), since a = 7/3. 

2.1. After SN Bursts Stop 

It is clear that an upper limit should be inserted into the basic equation (6); this is the 
time after which the bursting phenomena stops, t^^^^^ . 

Since there is the pdV term in the the first theorem of thermodynamics [dQ = = 
dU + pdV), the total thermal energy decreases with time. The pressure of the internal gas 
decreases according to the adiabatic law, 

P = Pi{^?'\ (11) 

where 

O T? j-burst 

= 3^^^ ' 

and the equation for the conservation of momentum becomes 

It is important to remember that this phase occurs after the SN burst stops. At that time 
{t = t*^™***), the volume of the bubble is computed by using i?j(t^"''***), and the expansion speed 
is equal to Rj{t^'^'^'^^), both of which are obtained from the numerical solution of equation (8). 
On assuming (also here) that p is constant irrespective of R or t, equation (12) is now inserted 
into equation (1) and the expansion law is obtained for the after-burst phase, 

fl.(iii)'"(^^s|!:!^y'\v. (14) 

which is identical to equation (10.33) of (McCray 1987). 

2.2. Numerical Integration 

The differential equations need to be solved by sectors, with each sector being treated as 
independent from the others, except for the coupling in computing the volume [ see equation (2)] 
of the bubble. From a practical point of view, the range of the polar angle 6 (180°) will be 
divided into uq steps, and the range of the azimuthal angle (360°) into steps. This will 
yield {ng +1) {n^ +1) directions of motion that can also be identified with the number of 
vertices of the polyhedron representing the volume occupied by the expansion; this polyhedron 
varies from a sphere to various morphological shapes based on the swept up material in each 
direction. 
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In 3D plots showing the expansion surface of the explosion, (see for example figure 12), 
the number of vertices is = {ne+l)-{n^ +1) and the number of faces is ne ■ n,j), typically we 
have 77,0=50 and n(/,=50 (in this case the subscript varies between 1 and 2601). However, all 
calculated models are axisymmetric, and the essential number of points to draw such figures is 
only 720 + 1=51. -Rup, -Req and -Rdown are now introduced, which represent the distances from 
the position of the OB associations (denoted by zqb) to the top, to the left and to the bottom 
of the bubble; see, for example, figure 2 where zq^ =100 pc. 

2.3. The Numerical Equations 

At each time step ,At, the volume V of the expanding bubble is computed [see equa- 
tion (2)]. In other words, the volume V swept up from the explosions is no longer a sphere, but 
becomes an egg or an hourglass. The pressure in the first phase, see equation (5), is computed 
through the following finite-difference approximation: 



/ = + 



7-1 pfe-i y'^ - 1/''"^ 



At, (15) 



V ' At 
where k is the number of steps considered. 
Equation (1) now leads to 

(p,i?j4) = 3pRl (16) 



dt 

which may rewritten as 



Si??/?/ + R^Mj = 3p^ - ^R^Ai. (17) 

' ' Pj Pj ' 

The first term represents the ram pressure of the stratified ISM on the expanding surface and 
the second represents the inertia of the bubble. The average density pj is numerically computed 
according to the algorithm outlined in subsection 2.5; the time derivative of pj at each time 
step At is computed according to the finite difference- approximation, 

(IB) 

where k is the number of steps considered. 

Equation (17) can be re-expressed in two differential equations (along each direction j) 
of the first order suitable to be integrated: 



dyi,j 
dt 



y2j, (19) 



—f = 3p- 3 ^ - ^|/2,i • (20) 

In the new regime {t > t^'^^^^^ equation (13) now becomes 
2 p V^^^ 

3R]Rj + R]Rj = 3 ^j./3^_ _ -^^R^.Rj. (21) 
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Equation (21) can be re-expressed as two differential equations (along each direction j) of the 
first order suitable to be integrated: 



dyi,j 
dt 



y2,j, (22) 



dy2,j _ SpiVi^^ _ _pj 

dt ~v'i'p,y., yi, p/"" ^ ^ 

The integrating scheme used is a Runge-Kutta method, and in particular the subroutine 
RK4 (Press et al. 1992); the time derivative of the density along a certain direction j, pj, is 
computed at each time step according to formula (18). 

The integration time, tage? and the time steps are always indicated in the captions 
connected with the various diagrams. 

From a numerical point of view, the pressure across the bursting time is continuous. 
This is because the first value of the pressure after the bursting time is that of the last step in 
the bursting phase, modified according to the adiabatic law modelled by formula (11). 

The upper limit chosen to integrate the differential equations is tage = 2.5 ■ 10^ yr; this is 
the approximate age of GSH 238, the super-shell simulated in subsection 4.2. 

2.4- The Density Profile 

The vertical density distribution of galactic HI is well-known; specifically, it has the 
following three component behaviour as a function of z , the distance from the galactic plane 
in pc: 

n{z) = riic-''^^^' + riae-"'/^^' + rige-l^l/^^ . (24) 

We took (Bisnovatyi-Kogan, Silich 1995; Dickey, Lockman 1990; Lockman 1984) ?t,i=0.395 
particles cm"'^, ifi=127 pc, n2=0.107 particles cm""^, i^^2=318 pc, ^3=0. 064 particles cm~^, 
and if3=403 pc. This distribution of galactic HI is valid in the range 0.4 < R < Rq, where Rq 
= 8.5 kpc and R is the distance from the galaxy center. A plot showing such a dependence of 
the ISM density from z is shown in figure 1. 

2. 5. Computation of the Swept Mass 

The ISM density is not constant, but varies in the ^-direction according to equation (24). 
The swept mass can therefore be computed in a certain direction j using the following algorithm: 

1. The pyramidal sector is divided into layers (for example 1000) whose radii range from 

Rl~1/2 to Rl+1/2- 

2. In each layer the volume AVl = |(-Ri+i/2 ~Rl-i/2) computed as well as the correspond- 
ing mass, AMl = AVLPj{z = Rl cos 9 j + zqb), where 9j represents an angle between the 
z-axis and the j-th radial path. 

3. The various contributions, AM^, are added in order to obtain the total swept mass in the 
considered sector. 
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Fig. 1. Average structure of the gaseous disk in tlie z -direction; z is allowed to vary between pc and 
400 pc. 

4. The average density along a sector j as pj is calculated using AM/, and AV^ as pj = 

2.6. Astrophysical Units 

Our basic units are: time (tj), which is expressed in 10^ yr units; £"51, the energy in 10^^ 
erg; rig the density expressed in particles cm~'^ (density po = nom, where m=1.4mH); and A^*, 
which is the number of SN explosions in 5.0 ■ 10^ yr. 

By using the previously defined units, formula (10) concerning "the expansion law" in 
the bursting phase becomes 



i? = 111.56 pc^^^ — )5. 

no 



(25) 



Conversely, equation (14) concerning "the expansion law" in the adiabatic phase is 



i? = 171.85 pci 



-^51^^i.,burst^ii 



no 



(26) 



This is the approximated radius derived by a spherical model. A perfect coincidence is not 
expected upon making a comparison with that expected in a stratified medium obtained with 
the multidimensional thin shell approximation. Another useful formula is the luminosity of the 
bubble (see McCray 1987), 

LsN = EoRs^ = 0.645=^^erg s-'E,,N*. (27) 

This formula can be useful in order to derive the parameter A^* based on observations; for 
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example, a luminosity of 1.6 lO^^erg s ^ corresponds to A^* = 248. The total deposited energy, 
^tot' is 

Etot = i^5iiV*10-^^erg, (28) 
when t''™''*=tagc and 

j-burst 

^tot = ^5iiV*10^^^-^erg, (29) 

The spectrum of the radiation emitted depends on the temperature behind the shock 
front [ see for example formula 9.14 in McKee (1987)], 

T-U^-l K». (30) 

where fi is the mean mass per particle, k the Boltzmann constant and vs the shock velocity 
expressed in cm sec~^. A formula which is useful for the implementation in the code is easily 
derived, 

T = 31.80 vl-^ K°, (31) 
when v^-^ is expressed in km sec~^. 
2. 7. Other Effects 

Two other effects are now analysed: the galactic rotation and the influence of the ambient 
pressure on the late phase of motion. 

2.7.1. Galactic Rotation 

The effects of galactic rotation are important in shaping the bubbles, as can be seen in 
figure 2 of Silich et al. (1996). This is very important on the galactic plane, and the properties 
of the distorted "walls" for a rotating galaxy were described by Palous et al. (1990). The 
application of this distortion to Gould's Belt was done by Moreno et al. (1999) , and the 
impact on large scales in a percolation-type scheme (in this case the cloud and star- formation 
mechanism is deterministic, not stochastic) was studied by Palous et al (1994). The distortion 
of the superbubble due to galactic rotation is now introduced. The rotation curve can be that 
given by Wouterloot et al. (1990), 

VniRo) = 220{^r''\msec-\ (32) 

where Ro is the radial distance from the center of the Galaxy expressed in pc. The translation 
of the previous formula to the astrophysical units adopted gives 
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^-f?o[pc]^ 0.382 



rad 



f](i?o) = 2244 ^ ^g^y ■ (34) 

Ro[pc\ W yr 

Here, fl{Ro) is the differential angular velocity and 

/ /?o[pc] \ 0.382 

^ = 2244 ^ 8500^ rad. (35) 

Here, 0, is the angle made on the circle and t-j the time expressed in 10^ yr units. Upon 
considering only a single object, an expression for the angle can be found once R — Rq -\- y is 
introduced, 

m = (36) 

Ro + y 

The shift in the angle due to differential rotation can now be introduced, 

A0 = 0(y)-0(O), (37) 

where x and y are the superbubble coordinates in the inertial frame of the explosion, denoted 
by x=0,y=0 and z=zob- The great distance from the center allows us to say that 
Ar 

^ = A0, (38) 

where Ax is the shift due to the differential rotation in the x coordinate. The shift can be 
found from (38) and (37) once a Taylor expansion is performed, 

Ax^-VniRo)^t. (39) 

Upon inserting (33) in (39), the following transformation, T^, due to the rotation is obtained 
for a single object in the solar surroundings: 

xf = x + 0.2Q4y t 

Tr I yf = y (40) 

zf = z, 

where y is expressed in pc and t in 10^ yr units. A typical distortion introduced on a circular 
section of the superbubble, z = zqb, is shown in figure 18. 

2.8. Late Evolutionary State 

The problem of the late evolutionary stages can be attached including the ambient 
pressure; but firstly the Euler number (ratio of the pressure forces to the inertial forces in a 
flow) should be introduced 

Eu = (41) 

where ps is the density behind the shock, Vs is the shock velocity and p the ambient pressure. 
The nominal value of pressure at the Solar circle, 10~^^dyne cm~^, can be used as an adequate 
value (see Boulares, Cox 1990). Another density distribution along the z-axis and a discussion 
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of the pressure distribution can be found in Kim et al. (2000). In order to simulate the influence 
of the ambient pressure, a drag acceleration should be introduced, 



dy2,j 



-Vd Eu, (42) 



dt 

where is Vd is an artificial coefficient. This effect is negligible when the Euler number is small 
(high shock velocities), and gives a negative contribution when the Euler number ^ 1. Only 
in subsection 4.1 are some results adopting the value of Vd = 6000pc/10'' yr ^ 600km sec~^ 
reported in comparison with the normal run, the associated note "with the Euler process" being 
always inserted. 

3. Test Calculation 

The shape of the superbubble depends strongly on the time elapsed since the first ex- 
plosion, tagc, the duration of SN burst , t^"''***, on the number of SN explosions in the bursting 
phase, and on the adopted density profile. 

In this section, where the initial tests are discussed the galactic rotation is neglected. 

Ordinarily, the dynamics of a superbubble is studied with hydro simulations (partial 
differential equations). Conversely, here under the hypothesis called "thin-shell approximation", 
the dynamics is studied by solving ordinary differential equations. This choice allows the 
calculation of a much larger number of models compared with hydro-dynamics calculations. 

3.1. Reliability of the code 

The level of confidence in our results can be given by a comparison with numerical hydro- 
dynamics calculations see, for example, MacLow et al. (1989). The vertical density distribution 
they adopted [ see equation(l) from MacLow et al. (1989) and equation(5) from Tomisaka, 
Ikeuchi (1986) ] has the following it z dependence, the distance from the galactic plane in pc: 

V{z) 



^hydro = rid { Oexp[-^^^^] + (1 - 6)exp 



o"ic 



J 



(43) 



with the gravitational potential as 
y(2) = 68.61n 



1 +0.9565 sinhM 0.758— 



(kms-^)^ (44) 



Here rid = 1 particles cm""^, O = 0.22, ajc = 14.4km s~^, ac = 7.1km s~^ and Zq = 124 pc. Table 1 
reports the results of ZEUS (see MacLow et al. 1989), a two-dimensional hydrodynamic code, 
when tagc = 0.45 ■ 10^ yr. The supernova luminosity is 1.6 ■ lO'^^erg s~^, zob=100 pc and the 
density distribution is given by formula (43). The ZEUS code was originally described by Stone, 
Norman (1992). 

In order to make a comparison our code was run with the parameters of the hydro-code 
(see figure 2); a density profile as given by formula (43) was adopted. 

From a practical point of view, the percentage of reliability of our code can also be 
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Fig. 2. Section of the superbubble on the x-z plane when the explosion starts at zqb = lOOpc. The density 
law is given by equation 43. The code parameters are iago=0-45 10^ yr, At ~ 0.001 • 10^ yr, tj^'^^^=0A5, and 
iV*=248. The points represented by the small crosses indicate the inner section from figure 3a of MacLow 
et al. (1989). 



where -Rhydro is the radius, as given by the hydro-dynamics, and Rnum the radius obtained from 
our simulation. In the aheady cited table 1, our numerical radii can also be found in the 
upward, downward and equatorial directions, and the efficiency as given by formula (45). The 
value of the radii are comparable in all of the three chosen directions, figure 2 also reports on 
the data from figure 3a of MacLow et al. (1989). 

Koo, McKee (1990) applied their one-dimensional method based on the virial theorem 
to a two-dimensional calculation by using a sector approximation, which is identical to that 
assumed in the present paper and in which a differential equation is solved along a radial path, 
figure 8 from (Koo, McKee 1990) is for the same density distribution as given by formula (43). 



introduced. 




hydro 



(45) 
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Table 1. Code reliability. 
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Fig. 3. Rcq/Rup versus the time, tage, expressed in units of lO'^ yr when the explosion starts at zqb = Opc. 
The parameters are At = 0.1 • lO'^ yr, t^"''''*= 0.45, and N*= 248. 

3.2. The Asymmetries 

Firstly, the time evolution of the ratio Rcq/Rup is focused on as a function of tage- 
figure 3 plots the asymmetry indicator, i?eq/-Rup for a model with zqb = 0, t^'^^^^ = 0.45 and 
A^* = 250. In this case, the indicator of the asymmetry, Rcq/R^p is ~ 0.9 at tage= 0.4- 10^ yr, 
reaches a minimum at tage= 2.0 ■ 10^ yr and remains nearly constant up to tage= 2.5 ■ 10^ yr (see 
figure 3 ). 

Secondly, we consider how the asymmetry indices, -Rdown/-Rup (figure 4) and Rcq/Rup 
(figure 6), depend upon the height of the OB association, zq^. A comparison is made at 
^age = 4.5- lO^yr and A^* = 248. Now, the first index of the asymmetry, -RdownZ-Rupj is 1 at 
2^06=0 pc and 0.52 at 2;ob=100 pc (value of zqb where the comparison with the hydro-code was 
performed; see figure 4). Another interesting theoretical parameter that has an observational 
counterpart is the degree of symmetry with respect to the galactic plane. Upon adopting the 
previously developed notation, two variables, Dup and -Ddown, are introduced when -Rdown > 2:03: 
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Fig. 4. i?down/i?up versus zqb in pc. The parameters are ta 
0.45, and N*^ 248. 



=0.45 • 10^ yr, At = 0.005 • 10^ yr, 



burst. 



-Run + ^OB — D 



up; 



-Rdown — ^OB —D 



down- 



(46) 



The super-shells are often visualised in galactic coordinates. When this visualisation is 
possible, the ratio between the positive maximum galactic latitude and the absolute value of 
the negative minimum galactic latitude belonging to a given super-shell can be identified with 
^""^ ■ In our theoretical framework, the ratio is easily found to be 

-'-^ c\ own 



D 



up 



up 



ZOB 



D R -z ^ 

-'-^down -'l-down ^OB 

and can be visualised in figure 5; due to the nature of the phenomena, this ratio easily takes 
high values. 

Conversely, the second index, i?eq/-Rup, as a function of zqb presents a marked minimum 
value at ^ 100 pc; for this value of zqb the asymmetry in the shape of the expanding envelope 
is maximum (see figure 6). 

3.3. Characteristic Shapes 

The developed theory of bubble expansion in the ISM has four key parameters: the 
distance zqb from the galactic plane , t^"''*, the bursting phase, A^*, the number of SN explosions 
in 2.5 • 10^ yr, and tage the age of the supershell . By varying these four parameters, a great 
variety of different shapes can be obtained and two cases, zob=0 pc and 2;ob=100 pc, will be 
concentrated on. 
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Fig. 7. Stars indicate tlie existence of "vertical walls"; iagc is expressed in 10'' yr units. The parameters 
arc A< = 0.1 • 10^ yr, zob=0, and t^"'''^* = iage- 

3.3.1. From the ellipse to the hourglass, zob=0 pc 

We start by fixing 2ob=0 pc (the symmetrical case in the +z and -z directions). If the 
luminosity is high enough and zq^ = 0, the superbubble changes its shape from spheroidal to 
hourglass shape via the vertical wall. In order to define the vertical wall, a reasonable angle 
of 30° is first introduced between the plane z = and a certain radius, -R300 ; this radius 
projected on the z=0 plane is denoted by -Rp,3o°. From this point of view, a vertical wall is 
made if Rz=o — Rp,m° Rz=o > Rp,30° indicates the elliptical shape and Rz=o < Rp,3o° indicates 
the hourglass shape. Another way to express the previous definition is the following inequality, 
which can be easily embedded in the code: 



-^eq R' 



p,30° I 



R, 



■ 100 < 5%. 



(48) 



■cq 



Through such a definition, the vertical walls will span 60° on the galactic height. The existence 
of the vertical wall can now be plotted in the tage — N* space; see figure 7 when t^urst^ ^^^^^ ^^^^ 
figure 8 when tj^'^^^=0.5. 

From these two plots (7 and 8) it can be inferred that: 

• When t^^^^^= 0.5, vertical walls are only developed if A^* > 20; for example at A^*=100 
they appear at tago ~ 1-0 ■ 10'^ yr (see figure 8). In this case the transition from an ellipse 
to an hourglass is not observed; in other words, the vertical walls are the typical structure 
over a wide range of tagc and A^*. 
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Fig. 8. Stars indicate the existence of the " vertical walls" ; tagc is expressed in lO^yr units. The parameters 
arc A< = 0.1 • 10^ yr, zoB=Opc, and i^"''"' = 0.5. 

• When t^i^fst^ f_^^^ transition between the various zones is better defined; for example, 
at A^*=50 the vertical walls appear at tagc ~ 1-0 • 10^ yr (see figure 7). 

A typical plot showing the temporal development of vertical walls through a displacement of 
the sections at regular time intervals is reported in figure 9. Further on, from the data of plots 
(7 and 8), we can extract an approximate age, after which the vertical walls appear when t^"''^*= 
0.5, and when t*^"''st= t^^gc, ^v^w- The numerical law can be obtained by fitting the value 
tvw (the value of tage at which the vertical walls appear) versus A^*, 

log(tvw) = a + &-log(iV*). (49) 

The values of a and b are easily derived from the least-squares fitting procedure. In other words, 
this relationship connects A^* and the first corresponding value of tage at the left of figure 7, 
and figure 8 is defined as tyw- The apphcation of least-squares gives 

(50) 
(51) 

It can therefore be concluded that introducing of a finite value of t^^"^^^ (for example t^"''***= 
0.5) inhibits the transition from the vertical walls to the hourglass shape. It is important to 
note that the time necessary to develop vertical walls decreases when the key parameter, A^*, 
increases; the exact law being an inverse power. In other words, the momentaneous radius of the 
evolving superbubble in the galactic plane (^=0) gives a reference radius; the balance between 
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Fig. 9. Section of the superbubble on the a;-z plane. The parameters are iago=2. 5- lO'' yr, Ai = 0.01 • 10^ yr, 
^OB=0pc. tj^'^^*'^ 0.5, and N*= 100. The sections at regular intervals of <ago/10 are also shown. 

a progressively greater distance at higher values of z and the decreasing matter encountered 
along the trajectory allows the radius along the other directions to be greater with respect 
to that along the galactic plane. Under the existence of certain parameters, as given by the 
numerical laws (50) and (51), a progressively increasing radius (connected with the progressive 
angle with respect to the galactic plane) having a "vertical shape" is obtained. 

3.3.2. Toward the V-shape, zqb > pc 

A substantial modification of the up/down symmetry with respect to the galactic plane 
can be obtained by varying the distance zqb from the galactic plane. In this case the mor- 
phologies become complicated and we limit ourselves to exploring the situation at 2;ob=100 pc 
where a transition from the egg shape to the V-shape is observed (see figure 10); in this case 
the temporal evolution of the superbubble is visualised through a displacement of the sections 
at regular time intervals. 

When 2;oB=100 pc the density in the equatorial plane of the explosions is ~ two-times 
lower with respect to the density of the galactic plane. The surface of the expansion will 
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500 



Fig. 10. Section of the superbubble on the x-z plane. The parameters are tago=2.5 • lO'' yr, 
At = 0.01 • 10'^ yr, zob=100 pc, t:^"''"*= 0.5, and N*= 250. The sections at regular intervals of tago/10 are 
also shown. 

therefore find an increasing density in the downward direction and, conversely, a decreasing 
density in the upward direction; the curious V-shape is generated from the competition of 
these two effects. 

The effect of increasing the bursting time is shown in figure 11, where tburst_ ^^^^ rather 
than = 0.5 ■ 10^ yr of figure 10. In the case of figure 11, the maximum elongation along 
the ^-direction is ~ 1400 pc against ~ 1100 pc of figure 10, in which a limited time of bursting 
was chosen. 

4. Single Astrophysical Objects 

Another way of detecting supershells is through a search for galactic worms; they are 
filamentary, vertical structures that look like worms crawling away from the galactic plane in 
HI and infrared surveys ( 60 and lOO/im) (Koo et al. 1992). 

With respect to the worms, there are also several alternatives to their origin. A possibility 
is due to gravo-magnetic instabilities, such as the Parker instability in 3D, which forms worm- 
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Fig. 11. Same as figure 10, but with i^urst^ ^^^^ 

like structures, as can be seen in figure 4b of Franco et al. 2002. 

Here, it is assumed tliat tliey are tlie vertical walls of supershells, as sustained by an 
observational point of view (Kim , Koo 2000) or by a theoretical argument (see appendix 1). 

The developed theory is now applied to a well-defined galactic supershell (see subsec- 
tion 4.1 and subsection 4.2 ). 

In this astrophysical section the galactic rotation is inserted. 

4-1- Supershell Associated with GW 46.4+5.5 

A careful study of the worms 46.4+5.5 and 39.7+5.7 (Kim , Koo 2000) has led to 
the conclusion that they belong to a single super-shell. Further on, the dynamical proper- 
ties of this HI supershell can be deduced by coupling the observations with theoretical ar- 
guments, Igumenshchev et al. (1990); the derived model parameters to fit observations are 
reported in table 2, where the altitudes of KK 99 3 and KK 99 4 (clouds that are CO emitters) 
have been identified with zqb by the author. 

These parameters are the input for our computer code (see the captions of figure 12). 
The problem of assigning a value to zqb now arises, and the following two equations are set up: 
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Table 2. Data of the supcrshcU associated with GW 46.4+5.5. 



Size (pc^) 


345 • 540 


Expansion velocity (km s"'^) 


15 


Age (10^ yr) 


0.5 


zos (pc) 


100 


Total energy (lO^^erg) 


15 



^up + zoB = 540 

1 ^° 

3. (52) 



Rnp + ZqB _ -Dup _ 15° 



-Rdown - ZOB -Ddown 3° 

The algebraic system (52) consists in two equations and three variables. The value chosen for 
the minimum and maximum latitude (15° and —5°) is in rough agreement with the position of 
the center at +5° of the galactic latitude (see Koo et al. 1992). One way of solving system 52 
is to set, for example, zob=100 pc. The other two variables are easily found to be as follows: 

zoB =100 pc, 

^down =235 pc, (53) 
-Rup =305 pc. 

For this value of zqb, we have the case where a transition from egg-shape to V-shape (as outlined 
in subsubsection 3.3.2) is going on, and the simulation gives the exact shape. In order to obtain 
Etot = 15. ■ lO^^erg (see formula (29)) and t*'"''''* = 0.5 ■ 10'^ yr, we have inserted A^*= 150. 

In order to test our simulation, an observational percentage of reliability is introduced 
that uses both the size and the shape, 

-Robs -Rnum 7 N /<- a\ 

eobs = 100(1 ' ^), (54) 

Z^j -'•'obsj 

where i?obs is the observed radius, as deduced by using the following algorithm. The radius 
at regular angles from the vertical (0°, 90°, 180°) is extracted from table 2, giving the series 
(305 pc, 172.5 pc, 235 pc); the cubic spline theory (Press et al. 1992) is then applied to 
compute the various radii at progressive angles (0°, 3.6°, 7. 2°,. ..180°), which is a series computed 
by adding regular steps of 180 ° /ng. The data are extracted from the dotted ellipse visible in 
figure 7 of Koo et al. 1992; this ellipse represents the super-shell at vlsr = 18.5 kmsec"^. 

We can now compute the efficiency over 50-1-1 directions [ formula (54)] of a section y-z 
when 2^0 which turns out to be eobs = 68.4% (eobs = 71.3% with the Euler process); the observed 
and numerical radii along the three typical directions are reported in table 3. The physical 
parameters adopted from (Kim , Koo 2000) turn out to be consistent with our numerical code. 
The result of the simulation can be visualised in figure 12, or through a section on the x-z plane 
(see figure 13). Another important observational parameter is the HI 21 cm line emission: the 
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Table 3. Radii concerning GW 46.4+5.5. 



Direction 


-Rnum(pc) 


^num(pc) with the Euler process 


-Robs(pc) 


Equatorial 


238 


233 


172.5 


Polar up 


342 


335 


305 


Polar down 


312 


304 


235 




Fig. 12. Model of GW 46.4+5.5. The parameters are iagc=0.5 • 10^ yr, At = 0.001 • 10^ yr, i^"^«t= 0.5, 
N*= 150, zoB=100 pc, and £'51 = 1. The three Eulerian angles characterising the point of view of the 
observer are $= 0° , 9= 90° and 0°. 

observation of H I gas associated with GW 46.4+5 (Hartmann, Burton 1997; Kim , Koo 2000) 
reveals that the super-shell has an expansion velocity of Kxp ~ 15 kms"^ (Kim , Koo 2000). 
In order to see how our model matches the observations, the instantaneous radial velocities 
are computed in each direction. A certain number of random points are then generated on the 
theoretical surface of expansion. The relative velocity of each point is computed by using the 
method of bilinear interpolation on the four grid points that surround the selected latitude and 
longitude (Press et al. 1992). 

Our model gives radial velocities of, Vtheo, 31 km s^^ < Vthco ^ 71 km s^^ (27 km s~^ 
< Vtheo < 66 km s~^ with the Euler process) and a map of the expansion velocity is given in 
figure 14, from which it is possible to visualise the differences in the expansion velocities among 
the various regions. 
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Fig. 13. Section of tlie supcrbubblc on tlic x-z plane when the physical parameters are the same as in 
figure 12. 

Perhaps it is useful to map the velocity (V^hco) the ^/-direction, 

^tLo = ^'(^,0)-sin(e)-(sin0). (55) 

This is the velocity measured along the line of sight when an observer stands in front of the 
super-bubble; 6 and are defined in subsection 2.2. The structure of the projected velocity, 
Vjheo' is mapped in figure 15 by using different colours; the range is km < V^^^ < 
36 km s~^ (0 km s~^ < V^hco — 32 km s~^ with the Euler process) and the averaged projected 
velocity is ~ 20km s~^ (~ 17.6km s^^ with the Euler process), a value that is greater by 
?a5km s~^ 2.6km s~^ with the Euler process) than the already mentioned observed expansion 
velocity, V^xp = 15 km s~^. 

As is evident from the map in figure 15, the projected velocity of expansion is not uniform 
over all of the shell's surface, but it is greater in the central region in comparison to the external 
region. In this particular case of an egg-shape, we observe a nearly circular region connected 
with the maximum velocities in the upper part of the shell. It is therefore possible to speak of 
egg-shape appearances in the Cartesian physical coordinates and spherical-appearances in the 
projected maximum velocity. 
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Fig. 14. Map of the expansion velocity relative to a simulation of GW 46.4+5.5 when 190000 random 
points arc selected on the surface. The physical parameters are the same as in figure 12 and the three 
Eulerian angles characterising the point of view of the observer are i>= ° , 9= 90 °, and \I'=0 °. 
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Fig. 15. Map of the velocity along the line of sight, V^^^, relative to a simulation of GW 46.4+5.5 when 
190000 random points are selected on the surface. The physical parameters are the same as in figure 15. 
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Fig. 16. Number of particles, iV^, along the line of sight relative to the simulation of GW 46.4+5.5 when 
NDIM=101. The other parameters are as in figure 15. 

Another way of presenting the results is through a ring-enhancement due to the law 
of e.m. emission in the external layer (see appendix 1). The projected number of particles 
N^{i,j), is then mapped using colour-coding contours (figure 16); the density enhancement in 
the ring region is evident from a visual inspection of this simulated map. 

Upon adopting formula (31), it is possible to map the field of temperature of the super- 
shell associated with GW 46.4+5.5 (see figure 17). 

The effect of galactic rotation on the shape of the superbubble is shown in figure 18 
where a circle is transformed into an ellipse. 

4.2. GSH 238 

The physical parameters concerning GSH 238, as deduced in Heiles (1998), are reported 
in table 4. The total energy is obtained based on a model derived by Heiles 1979 and the 
kinematic age is derived from approximate arguments. Once again, in order to obtain Etot = 
34. lO^^erg [ see formula (29) ] and t'^"''"* = 0.015 10^ yr, we have inserted A^*= 11300. 

The efficiency over 50+1 directions computed with formula (54) of a section x-y when 
z=0, is eobs = 58.4%. The 2D cut at z=0 of the superbubble can be visualised in figure 19. Our 
model gives a radial velocity at z =0 Vtheo=6km s~^ against the observed 8 km s~^ at 6=0. 
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Fig. 17. Map of the temperature relative to the simulation of the supershell associated with 
GW 46.4+5.5; the other parameters are the same as in figure 14. 

Table 4. Data of the super-shell associated with GSH 238. 

Size (pc^) 440 -1100 at b=0 

Expansion velocity (km s ^) 8 
Age (Wyr) 2.1 
Total energy (lO^^erg) 34 
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Fig. 18. Circular and rotation-distorted sections of the superbubble on the a;-?/ plane, denoted by z=zob- 
The physical parameters are the same as in figure 12; the rhombi represent the circular section and the 
stars the rotation-distorted section; the Galaxy direction of rotation is also shown. 

5. Collective Effects 

It should be remembered that the worms have an average linear extent in the z direction 
of 150 pc and a median width of 70 pc; the maximum extension that can be reached in the 
z-direction is ~ 1 kpc (Koo et al. 1992). These observations can be simulated by introducing 
a coupling between the percolation network reviewed in appendix 2 and the explosion model 
described in section 2. The percolation offers: 

1. The Cartesian coordinates of the explosions in the galactic plane that follow the grand- 
designed spiral galaxies. 

2. The times of the clusters vary between 10^ yr and 10^ yr. This cluster age divided by 4 can 
be identified with tage, the age of the superbubble. With this last choice, the maximum 
tage is 2.5 ■ 10^ yr, according to the age of GSH238, the oldest super-shell analysed here. 

In other words, the time sequence of the clusters (1,2, ,10-10^ yr) is converted in the 

following age sequence of the super-shells (0.25,0.5, ,2.5 -10^ yr). This choice fixes the 

maximum altitude reached by the super-shells in ^ 600 pc. 
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Fig. 19. Rhombi represent the circular section; the stars, the rotation-distorted section; and the small 
crosses, the observed data of GSH238 at b =0 (our z=0) as extracted from figure 8 of Heilcs (1998). The 
Galaxy direction of rotation is also shown. The parameters are tagG=2.1 • 10'' yr. At — 0.001 • 10^ yr, 



^bmst^ 0.015, N*= 11300, zoB=0 pc, and ^51 = 1- 



In order to test the suggested coupling between superbubbles and the percolation the distance of 
200 pc, which the explosions should cover in one time step, 10^ yr was inserted in equation (26), 
in order to trigger new cluster formations. 

Upon adopting the standard values of uq = 0.4particles cm~^, Eq=1 and A^*=150, we 

obtain t^urst^Q ggg^ ^g^l^g ^^g^^ 

seems to be in agreement with our previous assumption in 
which the bursting phase is lower than the age of the superbubble. A typical percolation run is 
shown in figure 20 for a face-on galaxy. The galactic rotation is considered in all the simulations 
discussed in this section. 

5.1. Few Superbubbles 

In order to visualise a few explosions of coordinates X and Y, only those clusters com- 
prised between XMIN[pc] <^ XMAX[pc] and YMIN[pc] <^ YMAX[pc] should be selected. The 
distance zqb from the galactic plane of the clusters is allowed to randomly vary between —100 
pc and 100 pc, which is the size of the active regions. A typical example of this network of 
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Fig. 20. Simulation of spiral galaxy Sh by adopting the theory and the parameters of appendix 2. The 
clusters are represented by spheres of decreasing radius with increasing age. The three Eulerian angles 
characterising the point of view of the observer are $= 25° , 8= 75°, and 25°. 



multi-explosions is shown in figure 21 through the algorithm that counts the number of particles 
along the line of sight; the selected boundaries are marked on plots and the appearance of the 
rings is evident. 

The superbubbles can also be visualised when the observer is situated in front of the 
galaxy plane (see figure 22). The effect of shape distortion due to formula (35) is evident. 

5.2. The Galactic Plane 

In order to simulate the structure of the galactic plane, the same basic parameters as in 
figure 20 can be chosen, but now the superbubbles are drawn on an equal-area Aitof projection. 
In particular, a certain number of clusters will be selected through a random process according 
to the following formula: 

selected clusters = pselect ■ number of clusters from percolation, (56) 

where pselect has a probability lower than one. The final result of the simulation is reported 
in figure 23. 
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Fig. 21. Number of particles, N^, along the line of sight relative to the network of the explosions when 
NDIM=200 and the galaxy is situated edge on. The boundaries of the box are XMIN=YMIN=-8000 pc, 
XMAX=YMAX=-4000 pc; 5 clusters were selected and each superbubble had 40000 random points on 
it's surface. The parameters were At = 0.01 • lO"^ yr, t^^rst^ q 2^ and N*= 50. 

We can also visualise the structure of the superbubbles as seen from an observer situated 
outside the Galaxy; the spiral structure arising from the ensemble of the shells is evident (see 
figure 24). The eUiptical shape of the superbubble, according to formula (35), is a function of 
both the age and the distance from the center. The distortion also follows the same inclination 
of the arms. 

6. Conclusions 

The expansion of a super-bubble in the ISM belonging to our galaxy can be simulated 
by applying Newton's second law to different pyramidal sectors. The following results were 
achieved: 

1. Characteristic structures, such as hourglass-shapes, vertical walls and V-shapes, could be 
obtained by varying the basic parameters, which are the bursting phase and the time over 
which the phenomena are followed. 

2. Single objects, like the super-shells associated with GW 46.4-1-5 and GSH238 were simu- 
lated with an efficiency eobs ~ 68% and eobs ~ 58%, respectively. 
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Fig. 22. Same as figure 21, but when the galaxy is situated face on. 




Fig. 23. Structure of the galactic plane in the Hammer-Aitof projection, as resulting from the super- 
bubble/percolation network. The value of pselect is 0.08, corresponding to 260 selected clusters. The 
parameters were At ~ 0.01 • 10^ yr, t7^'^^^=0.5, N*— 100 and each supcrbubblc had 200 random points 
on its surface. 
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Fig. 24. Number of particles, N^, along the line of sight relative to the network of the explosions when 
the galaxy is face on. Physical parameters were as in figure 23 and NDIM= 401. 

3. The network of many explosions that originate from the galactic plane could be tentatively 
simulated. 

Some problems are still unsolved: 

1. The results were obtained by using a standard law concerning the dependence of the ISM 
on the distance from the galactic plane, z. Perhaps a more refined treatment of the HI 
thin disk will change the obtained morphologies. 

2. Recent maps based on measurements of Nal absorption in the surroundings of the Local 
Bubble suggest the existence of narrow tunnels (Welsh et al. 2003). Can a tunnel be 
explained by the interaction of a super-bubble arising from the galactic plane (GSH 238) 
with a smaller interstellar bubble (the Local Bubble) ? 

Appendix 1. The Density Projection Algorithm 

Consider a shell within two concentric spheres with radii R and R + AR. The total in- 
tensity of thermal e.m. radiation, when absorption and scattering are neglected, is proportional 
to the following parameter (chosen among others): / oc / , where / is the dimension of the 
radiating region along the line of sight; /ring and /center denote, respectively, the intensity along 
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the ring and along the center. The length of the emitting region when the line of sight touches 
the inner shell is 



/ring = 2^(i?+Ai?)2-i?2. (Al) 

Conversely, when the line of sight is aligned with the center of the two shells the emitting length 
is 

^center = 2A_R. (A2) 

We briefly remember that on applying mass conservation to a region of radius R before and 
after an explosion and the junction condition concerning the density, the following formula 
concerning the depth of the expanding layer is obtained: 
AR 7-11 1 



R 7+13 12' 

where 7 = f (see Deeming, Bowers 1984). By using a little algebra, we find 



(A3) 



= 5. (A4) 

'center 

This argument is easily generalised in the case of AR = R/n^ where rit is an integer greater 
than 12 : 

Ir 



I, ring 



when A R = R/ut. (A5) 



center 



Therefore, the ratio /""^ = /""'^ can be called "the center darkening law", and says that the 

center ^center 

center of the superbubble should have a lower emergent intensity than the shell region (the 
limb). We point out that on the expanding surface of the superbubble there is a certain value 
of the intensity [for example, 30 in the band of 100 fiin (far-IR) (see figure 1 (gray-scale 
intensity) in (Kim , Koo 2000))]. The intensity at the center of the superbubble is theoretically 
predicted to be approximately five times lower than that of the ring (therefore, 5 ^^^). This 
argument was applied to our superbubble by plotting the density of random particles projected 
on the sky. In order to implement such a projection, we first computed the number of random 
elements {N{i,j,k)) that fall in a little cube of volume side^/NDIM^ where side is considered to 
be the length of the cube that encloses the superbubble, and NDIM the number of considered 
pixels. We then sum over one index, for example k, in order to obtain the elements N^{i,j) that 
fall in a given area — j^^pj^a • The already described algorithm can be considered a simulation 
of the number of particles that are being emitted on the line of sight. 

Appendix 2. Percolation and spiral galaxies 

The appearance of arms can be simulated through the percolation theory 
(Seiden, Gerola 1979; Seiden 1983; Schulman, Seiden 1986; Seiden , Schulman 1990); here, 
a previous model (Zaninetti 1988) will be slightly modified. The fundamental hypotheses and 
the parameters adopted in the simulation are now reviewed: 
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1. The motion of gas on the galactic plane has a constant rotational velocity, named Vq 
(in the case of spiral type Sb 218 kms~^). Here the velocity, Vq, is expressed in 200 
kms~^ units and will therefore be Vq =1.09. 

2. The polar simulation array made by rings and cells has a radius Rq = 12 kpc. The number 
of rings, (60), is the multiplication together of Rg and the number of rings for each kpc, 
named nring/kpc, which in our case is five . Every ring is then made up of many cells, 
each one with a size on the order of the galactic thickness, ~ 0.2 kpc. The parameter 
nring/kpc can also be found by dividing Ikpc by the approximate cell's size. 

3. The global number of cells, 11121, multiplied by the probability of spontaneous new cluster 
formation, psp (for example 0.01), allows the process to start (with the previous parameters, 
111 new clusters were generated) . Each one of these sources has six new surroundings that 
are labelled for each ring. 

4. In order to better simulate the decrease of the gas density along the radius, a stimulated 
probability of forming new clusters with a linear dependence by the radius ,Pst = a + bR, was 
chosen. The values a and b are found by fixing prmax ( for example 0.18), the stimulated 
probability at the outer ring, and prmin (for example 0.24) in the inner ring; of course, 
prminyprmax. This approach is surprisingly similar to the introduction of an anisotropic 
probability distribution in order to better simulate certain classes of spirals (Jungwiert , 
Palous 1994). 

5. Now, new sources are selected in each ring based on the hypothesis of different stimulated 
probabilities. A rotation curve is imposed so that the array rotates in the same manner 
as the galaxy . The procedure repeats itself n times (100); we denote this by Iq the age of 
the simulation, 100 ■ lO^yr being 10'' yr, the astrophysical counterpart of one time step. 

6. In order to prevent catastrophic growth, the process is stopped when the number of sur- 
roundings is greater than max (1000) and restarts by spontaneous probability. 

7. The final number of active cells (3258) is plotted with the size, which decreases linearly 
with the cluster age. In other words, the young cluster are bigger than the old ones. Only 
ten cluster ages are shown; only cells with an age of less than life (in our case 10 ■ 10'^ yr) 
are selected. A more deterministic approach has been followed by Palous et al (1994), 
where a model that connects the energy injection by star formation with the resulting 
interstellar structures in a differentially rotating disc has been introduced. 



References 

Basu, S., Johnstone, D., & Martin, P. G. 1999, ApJ, 516, 843 
Begelman, M. C. ,& Li, Z.-Y. 1992, ApJ, 397, 187 

Bisnovatyi-Kogan, G. S. ,& Silich, S. A. 1995, Rev. Mod. Phys., 67, 661 
Boulares, A. ,& Cox, D. P. 1990, ApJ, 365, 544 



35 



Deeming, T.,&: Bowers, R.L. 1984, Astrophysics II: Interstellar Matter and Galaxies, (Boston: Jones 

& Bartlett Pub) 
Dickey, J.M. ,& Lockman, F.J. 1990, ARA&A, 28, 215 
Ferriere, K. M., Mac Low, M., & Zweibel, E. G. 1991, ApJ, 375, 239. 
Franco, J., Kim, J., Alfaro, E. J., & Hong, S. S. 2002, ApJ, 570, 647 

Hartmann,D. ,& Burton, W.B. 1997, Atlas of Galactic Neutral Hydrogen (Cambridge: Cambridge 

University Press) 
Heiles, C. 1979, ApJ, 229, 533 
Heiles, C. 1998, ApJ, 498, 689 

Igumenshchev, L.V.,Sliustov, B.M. ,& Tutukov, A.V. 1990, A&A , 234, 396 
Jungwiert, B. ,& Palous,J. 1994, A&A, 287,55 
Kamaya, H. 1998, ApJL, 493, L95 

k Kim, J., Franco, J., Hong, S. S., Santillan, A., & Martos, M. A. 2000, ApJ, 531, 873 

Kim,K.T. ,& Koo,B. C. 2000, ApJ , 529, 229 

Koo, B. C. ,& McKee, C. F. 1990, ApJ, 354, 513 

Koo,B.C.,Heiles,C. ,& Reacli,W.T 1992, ApJ , 390 , 108 

Lockman, F.J. 1984, ApJ , 283, 90 

Mac Low, M.-M. ,& McCray, R. 1988, ApJ, 324, 776. 

MacLow, M.-M. , McCray, R. ,& Norman,M.L. 1989, ApJ , 337,141 

McCray, R. 1987, in Spectroscopy of Astrophysical Plasmas, ed. A. Dalgarno A. & D. Layzer 

(Cambridge: Cambridge University Press) 271 
McCray, R. ,& Kafatos , M. 1987, ApJ , 317,190 

McKee C.F., 1987, in Spectroscopy of Astrophysical Plasmas, ed. A. Dalgarno A. &: D. Layzer 

(Cambridge: Cambridge University Press) 226 
Moreno, E., Alfaro, E. J., & Franco, J. 1999, ApJ, 522, 276 
Palous, J., Franco, J., ,& Tenorio-Tagle, G. 1990, A&A, 227, 175 
Palous, J., Tenorio-Tagle, G., & Franco, J. 1994, MNRAS, 270, 75 
Pikel'ner, S. B. 1968, Astrophys. Lett., 2, 97 

Press, W. H. , Teukolsky, S. A. , Vetterling, W. T. , &; Flannery, B. P., 1992 Numerical Recipes in 

Fortran , second edition (Cambridge: Cambridge University Press) , 117 
Puche,D.,Westpfahl,D.,Brinks,E. ,& Roy,J.-R. 1992, AJ , 103, 1841 
Santillan, A., Franco, J., Martos, M., & Kim, J. 1999, ApJ, 515, 657 
Schulman, L. S. ,& Seiden, P. E. 1986, Science, 233, 425 
Seiden, P. E. 1983, ApJ, 266, 555 
Seiden, P. E. ,& Gerola, H. 1979, ApJ, 233, 56 
Seiden, P.E. ,& Schulman, L.S. 1990, Advances in Physics , 39,1 
Silich, S. A. 1992, Ap&SS, 195, 317 

Silich, S. A., Franco, J., Palous, J. ,& Tenorio-Tagle, G. 1996, ApJ, 468, 722 

Silich, S. A., Mashchenko, S.Y. , Tenorio-Tagle, G., ,& Franco, J. 1996, MNRAS, 280, 711 

Stone, J.M. ,& Norman, M. L. 1992, ApJS, 80, 753 

Tenorio-Tagle, G. ,& Bodenheimer, P. 1988, ARA&A, 26, 145 



36 



Tomisaka, K. 1992, PASJ, 44, 177. 
Tomisaka, K. 1998, MNRAS, 298, 797. 
Tomisaka,K. ,& Ikeuchi,S. 1986, PASJ, 38,697 

Walter, F., Kerp, J., Duric, N., Brinks, E. ,& Klein, U. 1998, ApJ, 502, L143 
Weaver, R., McCray, R., Castor, J., Shapiro, R, ,& Moore, R. 1977 , ApJ, 218, 377 
Welsh, B. Y., Lallement, R., Vergely, J. L., Crifo, F., & Sfeir, D. 2003, BAAS, 202,53 
Wouterloot, J. G. A., Brand, J., Burton, W. B., & Kwee, K. K. 1990, A&A, 230, 21 
Zaninetti, L. 1988, A&A , 190, 17 



37 



